Blantyre air quality modelling

Load required libraries

library(tidyverse)
library(here)
library(brms)
library(tidybayes)
library(sf)
library(mgcv)
library(priorsense)
library(osmdata)
library(nngeo) 

Air quality data

By summary, a random household survey was conducted in Blantyre Malawi during 2019-2020 (pre-COVID in Malawi) as part of a TB cluster randomised trial pre-intervention prevalence survey (sumamrised here: https://europepmc.org/article/MED/37862284). Fieldworkers carried purple air monitors attached to their backpacks, with measurements set to read every 90 seconds. Not all fieldworkers had a monitor. Previously, Helen Savage cleaned these monitor data, and merged to household questionnaire data, filtering measurements by interview times. Therefore, measurements are indoor household measurements.

Cluster boundaries (72 in total) are based on Ministry of Health Community Health Worker catchment areas, refined in collaboration with researchers. These were the unit of randomisation for the SCALE cluster randomised trial.

All air quality measurements were then taken from households within the 72 SCALE cluster.

load("input_data/aq_in.RData")
#input the scale clusters
load("input_data/scale_72_clusters.rda")

Get data into correct shape and aggregate measurements per household. Note as purple air devices took measurements every 90 seconds, mostly housholds surveyed have multiple measurements. Here, for convenience, we take the mean of these measurement per household.

Covariates

Population counts and density

#fix the coordinate crs of the scale clusters
scale_72_clusters <-  st_transform(scale_72_clusters, st_crs(mw_100m))

old-style crs object detected; please recreate object with a recent sf::st_crs()

Building density (i.e. we get the building footprints, and calculate their combined footprint area per grid cell, then calculate what percentage of the grid cell is “building footprint”) https://sites.research.google/gr/open-buildings/#open-buildings-download

Using Version 2 here - can update later to V3, with better precision, and year-specific data. (I think this is for 2020 currently)

#read in the data set
#is a stars dataset
buildings <- read_rds("input_data/blantyre_buff_buildings.rds")

#match the crs to the population data
buildings <-  st_transform(buildings, st_crs(mw_100m))

#convert stars object to polygons
mw_100m_grid_sf <- st_as_sf(mw_100m_cropped, as_points = FALSE, merge = FALSE)

#add grid cell ID for grouping
mw_100m_grid_sf <- mw_100m_grid_sf %>%
  mutate(grid_id = row_number())

#join building footprint data to the population data
buildings_with_grid <- st_join(buildings, mw_100m_grid_sf, left = FALSE)

#filter by confidence of building footprint, and by footprint size
#in Hannah R's previous experiements, this was a good compromise
#and matches well with home visits on the ground.
buildings_with_grid <- buildings_with_grid %>%
  filter(area_in_meters>20) %>%
  filter(area_in_meters<300) %>%
  filter(confidence>0.69)

#sum building footprint per grid cell
building_area_per_cell <- buildings_with_grid %>%
  group_by(grid_id) %>%
  summarise(total_building_area_m2 = sum(area_in_meters, na.rm = TRUE))

#cell area in m squared
cell_area_m2 <- 100 * 100

#compute percentage coverage
building_area_per_cell <- building_area_per_cell %>%
  mutate(building_coverage_pct = (total_building_area_m2 / cell_area_m2) * 100)

#join
mw_100m_grid_sf <- mw_100m_grid_sf %>%
  left_join(building_area_per_cell %>% st_drop_geometry(), by = "grid_id") %>%
  mutate(building_coverage_pct = replace_na(building_coverage_pct, 0))

#check by plotting
mw_100m_grid_sf %>%
  mutate(x = st_coordinates(st_centroid(.))[, 1],
         y = st_coordinates(st_centroid(.))[, 2]) %>%
  ggplot() +
  geom_tile(aes(x=x, y=y, fill=building_coverage_pct)) +
  scale_fill_viridis_c(option = "cividis") +
  theme_ggdist() +
  theme(panel.background = element_rect(colour = "grey78"))

Distance to the nearest main road. I.e. because we think that PM2.5 exposure can come from road vehicles.

We use Open Street Map (OSM) data to get the main roads (there are few of these in Blantyre… but see the most traffic). Later could consider including some smaller road…

#check plot
mw_100m_grid_sf %>%
  mutate(x = st_coordinates(st_centroid(.))[, 1],
         y = st_coordinates(st_centroid(.))[, 2]) %>%
  ggplot() +
  geom_tile(aes(x=x, y=y, fill=dist_to_road_m)) +
  scale_fill_viridis_c(option = "G", direction = -1) +
  theme(panel.background = element_rect(colour = "grey78")) +
  theme_ggdist() +
  theme(panel.background = element_rect(colour = "grey78"))

Warning: There were 2 warnings in `stopifnot()`.
The first warning was:
ℹ In argument: `x = st_coordinates(st_centroid(.))[, 1]`.
Caused by warning:
! st_centroid assumes attributes are constant over geometries
ℹ Run ]8;;ide:run:dplyr::last_dplyr_warnings()dplyr::last_dplyr_warnings()]8;; to see the 1 remaining warning.

Join all together into the modelling dataset, and just last tidy-up.

# Convert aq_data to sf object
aq_model_data <- st_as_sf(aq_model_data, coords = c("x", "y"), crs = st_crs(mw_100m_cropped))
aq_model_data <- st_transform(aq_model_data, st_crs(mw_100m_cropped))

aq_model_data <- st_join(aq_model_data, mw_100m_grid_sf)

aq_model_data <- aq_model_data %>%
  mutate(x = st_coordinates(st_centroid(.))[, 1],
         y = st_coordinates(st_centroid(.))[, 2]) %>%
  st_drop_geometry() %>%
  mutate(building_coverage_pct = building_coverage_pct/100)

aq_model_data <- aq_model_data %>%
  mutate(building_coverage_pct_z = scale(.$building_coverage_pct)[,1])

Model 1:

“We fitted a spatially smooth regression model for log-transformed PM2.5 concentrations using a Gaussian process (GP) smooth over spatial coordinates (x, y) to capture flexible spatial trends. Seasonal variation was modelled using a Fourier series expansion of day-of-year up to the 3rd harmonic (i.e., sine and cosine terms for annual, semi-annual, and tri-annual cycles) to account for complex seasonal patterns. The model was fitted in a Bayesian framework using brms with a Gaussian likelihood and the cmdstanr backend.”

Note, priors are still causing issues for me, so we will comment out now, and use the defauly stan weakly informative priors

Priors

priors <- c(
  prior(normal(3.4, 1), class = "Intercept"),         # Weakly informative on log_pm2_5
  prior(exponential(1), class = "sdgp"),              # GP standard deviation prior
  #prior(student_t(3, 0, 5.9), class = "sds"),               # Smooth term std dev prior (key!)
  prior(exponential(1), class = "sigma"),             # Observation error
  prior(normal(0,10), class = "b")
)

Fit prior only model

Now model with data.

Compare prior-only model to model with data using priorsense package:

# m1_draws <- as_draws_df(m1)
# 
# # Extract prior and posterior draws
# powerscale_sensitivity(m1, variable = c("b_Intercept", "b_sin_doy", "b_cos_doy", 
#                                         "sdgp_gpxy", "sigma", "Intercept"))
# 
# powerscale_plot_dens(m1, variable = c("b_Intercept", "b_sin_doy", "b_cos_doy", 
#                                         "sdgp_gpxy", "sigma", "Intercept"))

Predictions by space and time. Here, although we only have data within clusters, we will predict for cells outside of clusters based on covariates. We will also predict for the month for whcih we do not have data.

Now just draw 50 random coordinates, and predict by time to see whether we captured the trend.

Model 2: With covariates

Now we include grid level covariates of distance to the road, population density, and building footprint percent, along witht the mean temp and huidity on the day of measurement (from the purple air monitors)

Again, priors to be fixed later

Predictions by space and time

Sample coordinates, and plot over time

Compare models, using LOO CV

loo_compare(loo_m1, loo_m2)

   elpd_diff se_diff
m2    0.0       0.0 
m1 -115.8      18.6 

TODO

• PRIORS
• OTHER COVARIATES
• COMPARE MODELS WITH DIFFERENT BASIS FUNCTIONS FOR GP AND SPLINES
• EXCEEDANCES
• LINK TO MTB INFECTION PREVALENCE DATA, AND MODEL